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Abstract 

In this paper, we consider methods to compute the coefficients of in- 
terpolants relative to a basis of polynomials satisfying a three-term re- 
currence relation. Two new algorithms are presented: the first constructs 
the coefficients of the interpolation incrementally and can be used to up- 
date the coefficients whenever a nodes is added to or removed from the 
interpolation. The second algorithm, which constructs the interpolation 
coefficients by decomposing the Vandermonde-like matrix iteratively, can 
not be used to update or downdate an interpolation, yet is more nu- 
merically stable than the first algorithm and is more efficient when the 
coefficients of multiple interpolations are to be computed over the same 
set of nodes. 

1 Introduction 

In many applications, we are interested in computing the coefRcients of a poly- 
nomial interpolation of discrete data relative to a basis of polynomials pk (x) of 
increasing degree k = . . . n. In practical terms, this means that given n -\- 1 
function values fi at the n + I nodes Xi, i = . . .n, we want to compute the 
coefficients Ck, fc = . . . n of a polynomial gn{x) of degree n such that 

n 

9n{xi) = ^CkPkix) = fi, i = 0...n, (1) 

fe=0 

That is, the polynomial gn{x) interpolates the n + 1 function values fi at the 
nodes cc,. 
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If we are only interested in evaluating g„(a;) at different x, then the method of 
choice is Barycentric Lagrange Interpolation (Berrut and Trefethen 2004) , which 
avoids representing the interpolant in any specific base. Such coefficient-based 
representations are usefuU, however, if we are interested in computing other 
quantities such as the integral or derivative of the interpolant, its L2-iiorm 
and/or performing other operations on it such as transforming it to another 
interval. Aditionally, we may also be interested in updating the coefficients 
when new data is added or existing data is removed. In Gonnet (2010), such a 
representation is used in an adaptive quadrature routine for just these purposes. 

In the following, we will assume that the polynomials pi (x) of degree i can 
be constructed using a three-term recurrence relation, which we will write a^ 



(2) 



with 



Pq{x) = 1, p-i{x) = 0. 
Examples of such polynomials are the Legendre polynomials Pk {x) with 



ak 



2k -I 

or the Chcbyshcv polynomials Tk{x) with 



Pk = 0, 7fe 



"0 



1, Oik 



\. Pk 



0, 



k - 1 
2k- 1 



(3) 



(4) 



The coefficients Ci of Equation ([T]) can be computed solving the system of 
linear equations 



I po{xo) pi{xo) 
Po{xi) pi{xi) 



\ PaiXn) Pl{Xn) 

which can be written as 



• ■• Pn{xo) \ 

... Pn{xi) 
■■■ Pn{Xn) J 
p(")(.(") — f(") 



/ Co \ / fo \ 
Cl ./l 



(5) 



(6) 

The matrix p(") is a Vandermonde-like matrix and the system of equations 
can be solved in 0{tt') using Gaussian elimination. As with the computation of 
the monomial coefficients, the matrix may be ill-conditioned (Gautschi 1983). 

Bjorck and Pereyra (1970) present an algorithm to compute the monomial 
coefficients of an interpolation without the expensive and potentially unstable 
solution of a Vandermonde system using Gaussian elimination, by computing 
first the coefficients of a Newton interpolation and then converting these to 
monomial coefficients. This approach was later extended by Higham (1988) 



^ Gautschi (2004) and Higham (1988) use different recurrences which can both be trans- 
formed to the representation used herein. 
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to compute the coefRcients relative to any polynomial basis satisfying a three- 
term recurrence relation. Both algorithms compute the coefficients in O(n^) 
operations and in Higham (1990), both methods are shown to be numerically 
stable when a proppcr ordering of the nodes Xi, i — . . . n is used. 

In Section [21 we will re-formulate the algorithms of Bjorck and Pereyra 
and of Higham and extend them to update the coefficients after a downdate, 
i.e. the removal of a node, of an interpolation. In Section [3] wc present a new 
algorithm for the construction of interpolations of the type of Equation ([T]) based 
on a successive decomposition of the Vandermonde-like matrix in Equation ([5]). 
Finally, in Section 21 we will present some results regarding the efficiency and 
stability of both algorithms. 

2 A Modification of Bjorck and Pereyra's and 
of Higham's Algorithms Allowing Downdates 

Bjorck and Pereyra (1970) present an algorithm which exploits the recursive 
definition of the Newton polynomials 

TTkix) = (X - Xk-l)Trk-l{x). (7) 

They note that given the Newton interpolation coefficients a^, the interpo- 
lation polynomial can be constructed using Horner's scheme: 

qn{x) = an, qk{x) = {x - Xk)qk+iix) + ak, k = n -1...0 (8) 

where the interpolation polynomial is gn{x) = 90(2^)- 

They also note that given a monomial representation for qk{x), such as 

n—k 

qk{x) = Y,bf\' 

1=0 

then the polynomial qk-i{x) can be constructed, following the recursion in Equa- 
tion as 

qk-iix) = (x - Xk-i)qk{x) + ak-i 

n — k 

= {x - Xk-i) fcf ''a^' + Qfc-i 

i=0 

n— fc-t-1 n—k 

= E btW-Xk-iY.bl'^x^ + ak-i 

i=l i^O 
n — k 

= bi%x-'^+' + J2{bt\ - + -bi'^Xk-i + ak-i. (9) 

i=l 
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(fc-1). 



From Equation ^ we can then extract the new coefficients b] 

{ b[% i = n-k + l, 

bf'^^l C\-x._i6f\ l<^<n-k, (10) 
[ -b'^fl'^Xk-i + ak-i, i ^ 0. 

Higham (1988) uses the same approach, yet represents the Newton polyno- 
mials as a linear combination of polynomials satisfying a three-term recurrence 
relation. Using such a representation 



(x) = 5]cfV.(x) (11) 



he computes qk-i{x) by expanding the recursion in Equation ([8]) using the 
representation Equation (fTTj) : 



n — k 

Qk-lix) = {X - Xk-l) cf'^Pijx) + Ofc-i 

1=0 

n—k n—k 



= X] ''a;pi(a;) - Xk-i ^ cf^Piix) + ak-i 

1=0 i=0 

n—k n—k 

= X] '^i''^ {oLiPi+i{x) - Apj(x) 4- 7iPi_i(a;)) - Xk-i ^ c|''Vi(a;) + aj^l?) 

i=0 8=0 

Expanding Equation (jl2p for the individual Pk{x), and keeping in mind that 
P-_i(x) = 0, we obtain 

7i — k+l n—k 



qk-i{x) = ^ cf\ai-ipi{x) - ^ (xfe-i + I3i)pi{x) 

i=l 1=0 
n — — 1 

+ ^ c|^j:\7i+ipi(a;) -I- afc_i. (13) 



By shifting the sums in Equation ([T^ and re-grouping around the individual 
Pkix) we finally obtain 



(x) = C'^lk'^n-kPn-k+l{x) + [c-^\_^an^k-l ~ C^^lkixk-l + Pn-k)) Pn-k{x) 

J2 - {xk-, + 13,) + c|^i7«+i) Mx) 

i=l 

4'=' (xfc_i + /3o) + c'^i'hi + flfc-i ■ (14) 
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Higham then extracts the new coefficients c^^ from Equation ([T4)) as: 

cf}^ai-i, i = n - fc + 1, 

cf\ai-i - cf\xk-i + /30 + c*^\7j+i, I <i < n ~ k, 
-c\^'' {xk~i + /?o) + ''71 + O'k-i, i = 

In botli algorithms, the interpolating polynomial is constructed by first com- 
puting the divided differences 

a* = /N, ■ • ■ i^Q...n (16) 

and, starting with qn{;x) = a„, and hence Cg"'' = a„ or ^q"-* = a„, successively 
updating the coefficients per Equation or Equation ([TU)) respectively. 

Alternatively, we could use the same approach to compute the coefficients 
of the Newton polynomials themselves 

k 

T^kix) = ^77fV^(a;)• 
^=0 

Expanding the recurrence relation in Equation ([7]) analogously to Equa- 
tion ((TU), we get 

7rfe+i(a;) = lilf^ a},-pk+\{x) + {r]'^\a]^^x - vf\xk + h)^ Pk{x) 

i=l 

-~r,'i\xk+P^)+r,f'^li. (17) 
We initialize with r^Q*^^ = 1 and use 



(fc+i) 

■nl 



7?>_\a,_i, i = k + l, 

rif\ai-i - Tyf'^ [xk + A), « = fc, 

?7-^:\aj-i - vf^ (xk + ft) + Vi+ilt+i, l<i <k, 



(18) 



to compute the coefficients for ■Kk{x), k = \ ...n. Alongside this computa- 
tion, we can also compute the coefficients of a sequence of polynomials gu {x) of 
increasing degree k 

k 

9k{x) = ^cf\i{x) 

i=Q 

initializing with c^^ = oq, where the Oi are still the Newton coefficients com- 

(k) 

puted and used above. The subsequent coefficients q , fc = 1 . . . n are computed 
using 

i ~ 1 , (fc) n ^ ■ ^ ; ^ ' 
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This incremental construction of the coefficients, which is equivalent to ef- 
fecting the summation of the weighted Newton polynomials and is referred to 
by Bjorck and Pcreyra as the "progressive algorithm" , can be used to efficiently 
update an interpolation. If the coefficients ?7|"'' and c|"'' are stored and a new 
node Xn+i and function value fn+i sltc added to the data, a new coefficient 
a„+i can be computed per Equation (|16p. the coefficients ""^^ computed per 
Equation (fTSl) and, finally, the Cj-"'' updated per Equation ((T^ . resulting in the 
coefficients cf^^^^ for the updated interpolation polynomial gn+i{x). 



Algorithm 1 Incremental construction of g„ (x) 



vl)'^ < a-o - /3o, 7][^^ ^ ao (init 77(1)) 

for k — 1 . . .n do 

Vq <— 0, Vi ^ Xk {init v) 
for i = 2 . . .k do 

^ {{xk + /3i_i)wi_i - 7i_iu,_2) /ai-i {compute the Vi) 
end for 

Qk ^ v(0 : fc - lYc^'^^^'i {compute gk-i{xk)) 

TTfe v'''?/'^-' {compute TTk{xk)) 

ak ^ {fk - .9fe)/7rfc {compute ak, Equation H^)) 

c^'"') <— [c('^^^);0] + OfcTj^'^) {compute the new c^''\ Equation fig)) ) 

^(fe+i) ^ [0;a(0 : fc). *T7('=)] - [{xk +^(0 : fc)). * J?^^); 0] + 

[7(1 : k). * r7('')(l : fc);0;0] {compute the new r]'^''+^\ Equation p^]) 
13: end for 



We can re- write the recursion for the coefficients jy^-'^'' of the Newton polyno- 
mials in matrix- vector notation as 

where t('''+^) is the {k + 2) x (fc + 1) tri-diagonal matrix 



/ -/3o 71 



rp(fe + l) ^ 



"0 



72 



V 



afe-2 ~/3k-i Ik 
ak-i -h 

Oik ) 

and IgXA: is a (fc-f 2) x {k^V) matrix with Xk in the diagonal and zeros elsewhere 

/ Xk \ 

Xk 



V 



Xk 

/ 
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The vectors V'^) = (,7^ , , . . . , )T and = . . - A+i^V 

contain the coefScients of the fcth and (/c + l)st Newton polynomial respectively. 

Given the vector of coefficients c^") = {c^\c^^\ . . . , cl"'')^ of an interpola- 
tion polynomial gn{x) of degree n and the vector of coefficients of the 
{n + l)st Newton polynomial over the n + 1 nodes, we can update the inter- 
polation for a new node Xn+i and function value fn+i as follows: Instead of 
computing the new Newton interpolation coefficient a„+i using the divided dif- 
ferences as in Equation ()16p . we choose a„+i such that the new interpolation 
constraint 

gn+l{Xn+l) ~ g-niXn+l) + a-n+lT^n+liXn+l) = fn+1 

is satisfied, resulting in 

fn+l ~ gn{Xn+l) /on\ 
a„+l = -. r [^^J) 

which can be computed by evaluating (7„(x„+i) and 7r„+i(a;„+i). Note that since 
TTn+iixi) = for i = . . . n, the addition of any multiple of 7r„_|_i(a;) to gnix) 
does not affect the interpolation at the other nodes at all. This expression for 
a„+i is used instead of the divided difference since we have not explicitly stored 
the previous ai, « = . . . n, which are needed for the recursive computation of 
the latter. 

We then update the coefficients of the interpolating polynomial using 
and then the coefficients of the Newton polynomial using 

such that it is ready for further updates. Starting with ?7q*''' ~ 1 and n = 0, this 
update can be used to construct gn{x) by adding each Xi and /i, i = . . .n, 
successively. 

The complete algorithm doing just that is shown in Algorithm [TJ The ad- 
dition of each nth node requires 0{n) operations, resulting in a total of 0{n^) 
operations for the construction of an n-node interpolation. 

This is essentially the progressive algorithm of Bjorck and Pcrcyra, yet in- 
stead of storing the Newton coefficients a^, we store the coefficients ry,-"^"'^'' of the 
last Newton polynomial. This new representation offers no obvious advantage 
for the update, other than that it can be easily reversed: 

Given an interpolation over a set of 71 + 1 nodes Xi and function values fi, 
i = . . .n defined by the coefficients c,-"'' and given the coefficients rj^"''^^^ of 
the (n -|- l)st Newton polynomial over the same nodes, we will downdate the 
interpolation by removing the function value fj at the node xj. The resulting 
polynomial of degree n — 1 will still interpolate the remaining n nodes. 
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We start by removing the root Xj from the (n 
solving 



l)st Newton polynomial by 



^(n+l) 



for the vector of coefficients 77'^"^ Since Xj is a root of 7r„+i(cc), the system is 
over-determined yet has a unique solutioro. We can therefore remove the first 
row of (t("+^^ — lo^i) ^-iid the first entry of J7'-"+^\ resulting in the upper- 
tridiagonal system of linear equations 



\ 



72 



an -2 



-{xj + /3„_i) 

Otn-l 



In 

-{Xj + Pn) 
Oin 



/ (") \ 



(21) 



which can be conveniently solved in 0{n) using back-substitution. 

Once we have our downdated rj("\ and thus the downdated Newton poly- 
nomial 7r„(x), we can downdate the coefficients of gn{x) by computing 

5„_i(x) = gn{x) - a*7r„(x) 

where the Newton coefficient a* would need to be re-computed from the divided 
difference over all nodes except Xj. We can avoid this computation by noting 
that gn-i{x) has to be of degree n — 1 and therefore the highest coefficient of 



gn{x), ch^\ must disappear. This is the case when 



'In 







and therefore 



"■J ~ (n) ■ 
Vn 

Using this a*, we can the compute the coefficients of gn~i{x) as 



>-i) 



„(") 



(«) 

Cn (n) 
Vn 



1.. 



1. 



(22) 



The whole process is shown in Algorithm [51 The downdate of an 77-node 
interpolation requires 0{n) operations. 



■^Note that the nx (n + l) matrix (t("+^' — Iq^Jj) has rank n and the null space p(xj) = 
{poixj),pi(xj), . . . ,p„+i{xj))^ since for V = (t("+i) - IgXj^ p{xj), Vi = aiPij^x{xj)-(xj + 
0i)pi{xj) + "fiPi-i{xj) = by the definition in Equation ^ and the right-hand side 77'"+^' 



IS consistent. 
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Algorithm 2 Remove a function value fj at the node Xj from the interpolation 
given by the coefficients c*^"^ 

1: J?!"^ <— ??i'+i^V'^™ {compute rj^""' from rj^^+^'i using back- substitution) 

2. ^ + + I3n)vt^) /«„-! 

3: for i — n - 2. . .0 do 
5: end for 

6: Oj c^n^ /vii^^ [compute the coefficient Oj) 

7: c'^""^^ <— c'"-' — fljjy^"^ (compute the new coefficients c^"~^') 



3 A New Algorithm for the Construction of In- 
terpolations 

Returning to the representation in Equation ([S]) , we can try to solve the Vandermonde- 
like system of linear equations directly. The matrix has some special charac- 
teristics which we can exploit to achieve better performance and stability than 
when using Gaussian elimination or even the algorithms of Bjorck and Pereyra, 
Higham or the one described in the previous section. 

We start by de-composing the (n + 1) x (n + 1) Vandcrmonde-like matrix 
P(") as follows: 





/ 


\ 


pin) _ 


p(n-l) 


p(«) 




I q' 


Pn{Xn) J 



The sub-matrix P^" is a Vandermondc-like matrix analogous to P^"). The 
column p'"-* contains the ??th polynomial evaluated at the nodes Xi, i = Q . . . n— 1 



V Pn{Xn-l) J 

and the vector contains the values of the first ... n — 1 polynomials at the 
node Xn 

q""" = (Po(a;„),Pl(x„), . . . ,Pn-l{Xn)) ■ 
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Inserting this into the product in Equation ©, we obtain 



p(«-l) 


p(") 


q' 


Pn (*^n ) 



,(«-!) 



f(n-l) 



) / V C„ / \ fn / 



which, when effected, resuhs in the pair of equations 

p("-i)c("-i)+p(")c„ = f("-i) 



(23) 



where the vectors c^"~^^ = (cq, ci, . . . , c„_i)^ and f("~i) — (/q, /i, 
contain the first n coefficients or function values respectively. 

Before trying to solve Equation ()23p , we note that the columns of the matrix 
p("-i) contain the first ... n — 1 polynomials evaluated at the same n nodes 
each. Similarly, contains the same polynomials evaluated at the node a;„. 
Since the polynomials in the columns are of degree < n and they are evaluated 
at n points, p("~i) actually contains enough data to extrapolate the values of 
these polynomials at Xn ■ Using Lagrange interpolation we can write 



Tl-l 



9i 



3=0 



where the 



lf\x)^X\n-l) 



are the Lagrange polynomials over the first n nodes Xi, i = Q . . .n 
write Equation ([M]) as 

qT ^ £(n)p(n-l) 

where the entries of the 1 x ri vector are 



(24) 

(25) 

1. We can 
(26) 



The entries of can be computed recursively. Using the definition in 
Equation ()25p . we define 

YliXn -Xj). 
j=0 



w.„ 



and re- write 



as 



(27) 
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Using the previous tf^ and Wn-l^ we can re- write Equation (j27|) as 

("-1) 4"^'^ 



which, re-inserted into Equation (P?]). gives 



„(n — 1) — 1) 



Xji Xi Wfi~\ Xi Xyi — \ Xn Xi Wji^i 

for all i < n — 1. We then compute the last entry i = n — 1 using 

Jn) _ Wn 



Wn~l{Xn - Xn^i)' 



(28) 



(29) 



Therefore, starting with fg^-* = 1, we can construct all the 'E^''\ k = 2...n 
successively. Since, given £^''\ k = 1 . . .n the construction of each additional 
^("■-1-1) requires 0{n) operations, the construction of all the £^''\ k = l...n 
requires a total of 0{n^) operations. 

Returning to the Vandcrnionde-likc matrix, inserting Equation (j26p into 
Equation (|23p . we obtain 

p("-i)c("-i)+pWc„ = f("-i) 
^^"^P^"-')c("-i)+p„(a;„)c„ = /„. (30) 

Multiplying the first line in Equation (|30p with from the left and subtracting 
the bottom equation from the top one we obtain 

from which we can finally isolate the coefficient c„: 

£(")f(n-l) „ 



(31) 



Having computed c„, we can now re- insert it into the first Equation in Equa- 
tion ()30p . resulting in the new system 

p(n-l)j.(„-l) ^ _ p(„)^^ (32) 

in the remaining coefficients c^""^^. 

Applying this computation recursively to P^"^ p("-i)^ . . . ^ p(i) we can com- 
pute the interpolation coefficients c'-"-' . The final coefficient cq can be computed 
as 

The complete algorithm is shown in Algorithm [3] Since the construction of the 
fc = 1 . . .71 requires 0{n'^) operations (for each £^^\ 0{k) operations are 
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Algorithm 3 Direct construction of gn{x) 



for i = 2 . . . n do 

pW ^ ((x + /3,_i).*p(''-i) -7,_ip(*-2))/a,_i (/?Z/P(")) 
end for 

£q"'"'' <— 1, wi <— xi — xo (mif ^q"""-* and wi) 

for i = 2 . . . n do 

<— 1 (construct wt) 

for j . . . i — 1 do 

^ - Xj) 

end for 

for j = . . . i — 2 do 

i Xi—Xj Wi^i 

end for 

< T 

end for 

for i = n . . .1 do 



(compute the , Equation \28\) ) 
(compute £^j[^i, Equation i29\) ) 



a i~ i(i)pil-i)°o-.t-i)-pi{xi) (compute coefficient c,, Equation ^) 

18: f -i— f — CiP*-*-* (update the right-hand side f) 
19: end for 

20: Co <— /o/p'"''(0) (compute the final Cq) 



required to compute Wk and 0(k) are required to compute the new entries) and 
the evaluation of Equation ()3ip requires 0(k) operations for each Ck, k ~ n . . .0, 
the total cost of the algorithm is in 0(n'^) operations. 

Note that, as opposed to the algorithm presented in Section [21 this algorithm 
can not be extended to update or downdate an interpolation. It has an advan- 
tage, however, when multiple right-hand sides, i.e. interpolations over the same 
set of nodes, are to be computed. In such a case, the vectors f'-''-', k = 1 . . .n 
need to be computed only once (Lines [S] to [12] of Algorithm [3]). For any new 
vector f, only the Lines [TBI to need to be re-evaluated. 



4 Results 



To assess the stability of the two new interpolation routines described herein, 
we will follow the methodology used by Higham (1988). Higham defines a set 
of interpolations consisting of all combinations of the nodes 



Al: 




A2: 


Xi 


A3: 




A4: 





~ cos(i7r/n), 

-cos[(i-Hi)7r/(n-|-l)] , 

-1 + 2i/n, 

ijn, 



(extrema of T„(a;)) 
(zeros of r„+i(a;)) 
(equidistant on [—1,1] 
(equidistant on [0, 1]) 
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with the right-hand sides 



Fl: 




F2: 


/ 


F3: 





(1,0,...,0)T 
1/(1 + 25.t2). 

for I = ... n. 

To avoid instabihties due to unfortunate orderings of the nodes Xi , the nodes 
and corresponding function values were re-ordered according to the same per- 
mutation that would be produced by Gaussian elimination with partial pivoting 
applied to the Vandcrmondc-likc matrix, as described in (Higham 1990). This 
ordering is optimal for the Bjorck-Pereyra and Higham algorithms and produces 
good results for the two new algorithms described herein. 

For each combination of nodes and right-hand sides, we compute, following 
Higham, the coefficients c for the Chebyshev base (see Equation ^) for n = 5, 
10, 20 and 30 and compute the quantities 



ERR — — - — - — , RES 



If -Pell 



where c* is the exact solution and u is the unit roundofl^ as defined by Golub 
and Van Loan (1996, Section 2.4.2). Note that for the special case of Al or A2 
using the Chebyshev base, the coefficients can also be computed efficiently and 
reliably using the Fast Fourier Transform (Battles and Trefcthcn 2004). 

Results were computed using Gaussian elimination (Gl|f|), Higham's exten- 
sion of the algorithm of Bjorck and Pereyra (BP/b(^), the incremental Algo- 
rithm [1] (INCR) and the direct Algorithm [3] (DIRECT). The exact values were 
computed in Maple (Char, Geddes and Gonnet 1983) with 50 decimal digits of 
precision using the interp function therein. 

Results were also computed for the interpolation downdate (DEL) described 
in Algorithm [2j Starting from c* and rj* , the exact coefficients for the interpo- 
lation gn{x) and the Newton polynomial T:n+i{x) respectively, we compute the 
coefficients c^""^^ and rj'"' for gn-i{x) and 7r„(a;), resulting from the removal 
of the rightmost function value fk at Xk, k = argmax^ Xi. The exact coefficients 
c* after deletion were computed and used to compute the quantities ERR and 
RES. 

The results are shown in Tables [11 to [T^ For each n, the largest values for 
ERR and RES are highlighted. For the problem sets over the nodes Al and A2 
(Tables [1] to [5]) , the condition of the Vandermonde-like matrix is always < 2 
(Gautschi 1983), resulting in very small errors for Gaussian elimination. The 
Bjorck-Pereyra/Higham algorithm generates slightly larger residuals than both 
the incremental and direct algorithms for both sets of nodes. The values for 



^AU results were computed using IEEE 754 double-precision arithmetic and hence u Ri 
2.2 X 10-1^. 

''For the tests in this section, Matlab's backslash-operator, which uses partial pivoting, was 
used. In cases where the matrix is rank-deficient, a minimum-norm solution is returned. 
^Algorithm 1 in (Higham 1988) was implemented in Matlab. 
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GE BP/H INCR DIRECT DEL 



n 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


0.00 


0.50 


2.00 


3.94 


3.20 


4.30 


0.00 


0.50 


0.00 


0.72 


10 


3.20 


5.32 


1.20el 


2.79el 


7.76 


2.21el 


2.26 


7.00 


0.00 


2.47 


20 


7.28 


2.16el 


1.61e2 


5.27e2 


8.92 


3.64el 


9.14 


2.60el 


0.00 


3.40el 


30 


2.61 


1.03el 


6.72e2 


2.65c3 


2.08el 


l.lle2 


3.51 


i.eoei 


0.00 


8.76el 



Table 1: Results for problem Al/Fl. 



GE BP/H INCR DIRECT DEL 



n 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


0.97 


2.40 


1.82 


3.19 


0.93 


2.07 


0.88 


1.68 


0.00 


1.70 


10 


2.22 


4.19 


1.37cl 


3.77cl 


4.94 


1.32el 


1.93 


3.00 


0.00 


3.80 


20 


2.11el 


8.96 


9.93cl 


3.47e2 


2.24el 


4.80el 


1.80el 


1.27el 


0.00 


5.39el 


30 


3.63el 


1.36el 


1.27e2 


4.84e2 


5.55el 


2.27e2 


4.31el 


2.88el 


0.00 


1.19e2 



Table 2: Results for problem A1/F2. 
GE BP/H INCR DIRECT DEL 



n 


ERR RES 


ERR RES 


ERR RES 


ERR RES 


ERR RES 


5 

10 
20 
30 


1.26 2.67 

2.12 5.30 

1.13 4.81 
1.99 7.35 


1.16 2.27 

7.27 1.66el 

8.28 2.68el 
6.33 2.53el 


1.26 2.70 
1.83 3.63 
1.78 6.74 
1.14 5.86 


1.43 2.99 

1.19 1.87 
3.02 l.Olel 
2.55 1.04el 


0.00 1.18 
0.32 2.67 
0.05 5.25 
0.61 6.09 



Table 3: Results for problem A1/F3. 



GE BP/H INCR DIRECT DEL 



n 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


3.55 


0.70 


5.07 


7.31 


3.79 


1.58 


4.96 


4.96 


1.48 


3.16 


10 


1.19el 


3.65 


2.03el 


3.31cl 


8.34 


1.08el 


1.16el 


2.24 


7.23 


7.47 


20 


1.66el 


1.23el 


4.64el 


1.63e2 


5.57el 


1.41e2 


1.61el 


2.58 


1.69el 


2.41el 


30 


6.48el 


2.28el 


1.24e2 


4.36e2 


4.45el 


2.29e2 


6.49el 


5.98 


4.71el 


9.82el 



Table 4: Results for problem A2/F1. 



GE BP/H INCR DIRECT DEL 



n 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


2.57 


3.29 


2.88 


4.52 


2.35 


0.97 


2.30 


1.20 


0.60 


3.25 


10 


5.27 


3.89 


1.67cl 


4.59cl 


4.94 


5.23 


3.94 


3.75 


4.65 


1.03el 


20 


8.40 


8.80 


5.09cl 


1.63e2 


3.89el 


9.85el 


8.44 


4.04 


1.15el 


2.43el 


30 


3.39el 


1.95el 


1.17e2 


4.31e2 


3.00el 


1.99el 


3.42el 


2.55el 


3.27el 


1.13e2 



Table 5: Results for problem A2/F2. 
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GE 


BP/H 


INCR 


DIRECT 


DEL 


n 


ERR 


RES 


ERR RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


1.44 


0.11 


1.40 1.71 


1.17 


1.45 


1.12 


1.67 


0.00 


2.76 


10 


2.73 


3.22 


6.06 1.14cl 


2.86 


3.75 


3.64 


6.34 


6.56 


7.47 


20 


1.52 


5.59 


7.34 2.42cl 


2.06 


6.34 


3.92 


1.34el 


1.27el 


1.81el 


30 


2.81 


1.19el 


7.24 3.08cl 


1.65 


6.07 


2.79 


9.46 


1.41el 


2.94el 



Table 6: Results for problem A2/F3. 
GE BP/H INCR DIRECT DEL 



n 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


1.41 


0.58 


6.12 


l.Olel 


2.04 


0.97 


1.48 


0.73 


0.60 


0.78 


10 


2.99 


2.50 


l.lSel 


2.00el 


2.16 


1.52 


2.62 


3.11 


0.50 


1.42 


20 


2.26e3 


4.85 


2.01el 


7.19el 


3.58el 


1.02el 


6.23el 


3.97 


0.55 


2.77 


30 


1.39c6 


9.51 


3.90el 


1.34e2 


5.42e4 


2.60el 


3.07e2 


7.33 


0.55 


2.20 



Table 7: Results for problem A3/F1. 
GE BP/H INCR DIRECT DEL 



n 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


1.02 


1.98 


2.19 


3.53 


0.69 


0.80 


0.73 


1.40 


0.40 


0.86 


10 


4.91 


7.58 


3.83 


1.03cl 


1.05 


1.82 


1.61 


2.46 


0.47 


0.85 


20 


3.65e3 


1.07el 


8.32 


2.54cl 


5.73e2 


2.94 


1.29 


4.25 


0.63 


2.26 


30 


4.02e5 


9.93 


3.23el 


1.09c2 


1.35e5 


1.25el 


4.98 


7.32 


0.65 


3.48 



Table 8: Results for problem A3/F2. 
GE BP/H INCR DIRECT DEL 



n 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


1.48 


0.86 


1.87 


1.98 


1.78 


1.89 


1.38 


0.36 


0.00 


2.23 


10 


2.31 


3.28 


8.96 


2.76el 


2.00 


2.97 


3.27 


3.80 


0.45 


1.12 


20 


2.30e3 


5.53 


3.81el 


7.38el 


1.12e2 


1.24el 


3.13el 


5.11 


0.47 


3.39 


30 


1.39e6 


l.lOel 


2.28e2 


1.88e2 


5.29e4 


2.69el 


2.41e2 


4.78 


0.49 


2.29 



Table 9: Results for problem A3/F3. 





GE 




BP/H 


INCR 


DIRECT 


DEL 


n 


ERR 


RES 


ERR RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


2.01e2 


0.55 


0.55 0.76 


5.40c2 


0.74 


1.57el 


0.48 


3.06 


3.80 


10 


2.49e4 


0.84 


0.45 0.76 


3.56e6 


1.67 


8.20e2 


0.60 


1.94 


3.55 


20 


2.77el5 


2.53 


1.87 1.69 


5.08el4 


1.78 


6.00e7 


1.63 


5.93 


1.41el 


30 


4.50el5 


0.00 


4.44 1.16 






3.47ell 


1.60 


8.93 


3.82el 



Table 10: Results for problem A4/F1. 
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GE 




BP/H 


INCR 


DIRECT 


DEL 


n 


ERR 


RES 


ERR RES 


ERR RES 


ERR 


RES 


ERR 


RES 


5 


1.82e2 


0.48 


0.55 0.65 


3.64c2 0.37 


1.12el 


0.45 


1.34 


2.27 


10 


1.35e5 


0.63 


0.40 0.85 


3.77e6 1.55 


6.28e2 


0.75 


0.99 


2.01 


20 


2.81el5 


3.64 


0.71 1.20 


6.57el4 0.93 


4.05e7 


1.57 


1.83 


3.30 


30 


4.50el5 


0.00 


0.40 1.27 




1.37ell 


1.73 


2.88 


1.27el 



Table 11: Results for problem A4/F2. 





GE 




BP/H 


INCR 


DIRECT 


DEL 


n 


ERR 


RES 


ERR RES 


ERR 


RES 


ERR 


RES 


ERR 


RES 


5 


2.41e2 


0.23 


3.96el 0.72 


8.24e2 


0.73 


8.23 


0.62 


9.95 


1.13el 


10 


3.79e5 


1.38 


1.55e3 0.88 


4.94e6 


1.35 


3.04e2 


0.55 


0.26 


0.28 


20 


2.81el5 


3.74 


4.84e6 1.21 


9.48el4 


0.87 


3.56e7 


1.03 


1.45 


2.58 


30 


4.50el5 


0.00 


1.02ell 1.47 






l.lSell 


1.49 


4.56 


1.97el 



Table 12: Results for problem A4/F3. 



ERR, however, are usually within the same order of magnitude for the three 
algorithms. 

For the nodes A3 (Tables [7] to [9|), the condition number of the Vandermonde- 
like matrix is 5.11e6 for n = 30, resulting in the errors of approximately that 
magnitude when Gaussian elimination is used. In general, both the Bjorck- 
Pereyra/Higham and the direct algorithm generate smaller errors and residues 
than Gaussian elimination. The errors for the incremental algorithm are due to 
cancellation while evaluating gn{xn+i) for Equation (|20l) since the intermediate 
coefficients c*^'^) are several orders of magnitude larger than the resull[l. 

Finally, the condition number of the Vandermonde-like matrix for the nodes 
A4 is 4.26el6 for n = 30, making it numerically singular and thus resulting in the 
complete failure of Gaussian elimination. Note that since in such cases Matlab's 
backslash-operator computes the minimum norm solution, the resulting residual 
error RES is quite small. For the first two right-hand sides Fl and F2, the 
Bjorck-Pereyra/Higham algorithm performs significantly better than the two 
new algorithms, since the magnitude of the intermediate coefficients does not 
vary significantly. For the right-hand side F3, however, the errors are larger, 
caused by truncation in computing the Newton coefficients a^. The incremental 
algorithm fails completely for all right-hand sides since the intermediate and 
final coefficients c''^'-' , k < n, are more than ten orders of magnitude larger than 
the function valuet0 and the numerical condition of gn{xn+i) in Equation (j20p 

^ In (Higham 1988), Higham shows that the coefficients can be written as the weighted 
sum of any of the intermediate coefficients c'"' = JUj A'j'^j*' i where the fij depend only on 
the nodes and the coefficients of the three-term recurrence relation. If the fij are 0(1) and 

the intermediate c^'"' are much larger than the c'"\ then cancellation is likely to occur in the 
above sum. 

lie* II = 2.23el3 for F3 and n = 30. 
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exceeds machine precision, resulting in numerical overflow. These relatively 
large coefficients also cause problems for the direct algorithm when evaluating 
the right-hand side of Equation ([5^ . where the original function values are 
clobbered by the subtraction of the much larger p(")c„. 

The errors and residuals for the downdate algorithm are shown in the right- 
most columns of Tables [T] to [T^ In general the errors of the downdate are 
relatively small for all test cases. The larger residues, e.g. for A2/F2, arc due 
to cancellation in the final subtraction in Algorithm [51 Line [T] 

5 Conclusions 

We have presented here two new algorithms for the construction of polyno- 
mial interpolations. The first algorithm (Algorithm [1]) offers no substantial 
improvement over that of Bjorck-Pereyra/Higham except that it can be eas- 
ily downdated. The second algorithm, which does not allow for updates nor 
downdates, is slightly more stable than the other algorithms tested and is more 
efficient when multiple right-hand sides need to be computed over the same set 
of nodes. 
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